library(cem)
library(tidyverse)

rm(list=ls())
load("~/data/virginclip_pop.RData")

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW==0)

# outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
#virginclip_pop$did2005 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2005

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
#virginclip_pop$did2004 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$did2018 <- virginclip_pop$forestha_2018 - virginclip_pop$forestha_2011
#virginclip_pop$did2003 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2003

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2018 <- virginclip_pop$did2018 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre

##### Coarsened Exact Matching: iterations #####

##### Model 1 #####

vars1 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
          "slope", "elev")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars1 = setdiff(colnames(virginclip_pop), vars1)
mat1 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars1, keep.all=TRUE)

##### Model 2 #####

vars2 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
          "slope", "elev", "distcities")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars2 = setdiff(colnames(virginclip_pop), vars2)
mat2 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars2, keep.all=TRUE)

##### Model 3 #####

vars3 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
          "slope", "elev", "pop_2005", "pop_2010")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars3 = setdiff(colnames(virginclip_pop), vars3)
mat3 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars3, keep.all=TRUE)


##### Model 4 #####

vars4 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
          "slope", "elev", "distroads")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars4 = setdiff(colnames(virginclip_pop), vars4)

# relax distance cuts
distroads_cut = seq(0,300,25)

mat4 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars4, 
             cutpoints = list(distroads = distroads_cut), keep.all=TRUE)


##### ESTIMATION STEP #####

est1 = att(mat1, ddd2017 ~ moratorium + pop_2000 + pop_2005 + pop_2010 + 
             agc_bcc_mean + distroads, data = virginclip_pop)

est2 = att(mat11, ddd2017 ~ moratorium + dist_palms_km + dist_timber_km + dist_logging_km + distroads +
             agc_bcc_mean + pop_2000 + pop_2005 + pop_2010, data = virginclip_pop)

est3 = att(mat17, ddd2017 ~ moratorium + dist_palms_km +  dist_timber_km + dist_logging_km + distroads +
             agc_bcc_mean, data = virginclip_pop)

est4 <- att(mat4, ddd2017 ~ moratorium, data = virginclip_pop)


# remove data and save results

rm(virginclip_pop)
save.image("~/results/cem-results.RData")
